function DecodingProbability_v9()

% clear workspace
clear all
close all
clc

% profile on
tic

T=40; % transmitted packets
T_min=1;
T_step=1;
q=2; % field cardinality
g=[0.5 0.5]; %  layer selection probability
k1=10;
k2=10;
K1=k1;
K2=k1+k2;

% Decoding 1. Layer(working?)
l1_prob=zeros(T,1);
for pkts_recv = T_min:T_step:T % for each nb of transmitted packets
    
    l1_prob(pkts_recv)=MatrixFun(pkts_recv,k1,k1,g(1),q);
    
%     disp(['Layer 1: ' num2str(pkts_recv) ' out of ' num2str(T)])
end

% Decoding 2. layer
l2_prob=zeros(T,1);
for pkts_recv = T_min:T_step:T % for each nb of transmitted packets
% for pkts_recv = 100 % for each nb of transmitted packets
    
    for i=0:K1
        if pkts_recv<K2
            l2_prob(pkts_recv)=0;
        else
            val1=MatrixFun(pkts_recv,K1,i,g(1),q);
            val2=MatrixFun(pkts_recv,K2-i,K2-i,g(2),q);
            
%             [pkts_recv i val1 val2]
            
%             l2_prob(pkts_recv)=l2_prob(pkts_recv)+val1*val2;
            l2_prob(pkts_recv)=l2_prob(pkts_recv)+val2;
        end
        
        % Mean the shit
%         l2_prob(pkts_recv)=l2_prob(pkts_recv)/(K1/2);
        
    end
    
%     disp(['Layer 2: ' num2str(pkts_recv) ' out of ' num2str(T)])
    
end

%
l2_prob

% Replace 0 with NaN
for k=1:length(l1_prob)
    if l1_prob(k)==0
        l1_prob(k)=NaN;
    end
    
    if l2_prob(k)==0
        l2_prob(k)=NaN;
    end
        
end

% Plotting
figure(1)
hold('on')
plot(1:T,l1_prob,'-*')
plot(1:T,l2_prob,'-*')
hold('off')
grid('on')
pbaspect([2.5 1 1])
set(gca,'XTick',0:T/10:T)
set(gca,'YTick',0:0.1:1)
xlim([0 T])
ylim([0 1])

% Save plot
% print(gcf,'uep_ew_analytic.eps')

% profile off

toc

end


function Pr = MatrixFun(pkts_recv,layer_length,rank,g,q)

Pr=0;

if pkts_recv>=rank % if we have less pkts_recv than required rank, no need to calculate -> impossible
    
    bino_vector=binopdf(1:pkts_recv,pkts_recv,g); % binomial vector with the probability of 1,2,3,...,pkts_recv from layer
    
    for outcome=1:length(bino_vector) % for each outcome of the received packets
        if outcome>=rank
            Pr=Pr+bino_vector(outcome)*ProbMatricesWithRank(outcome,layer_length,rank,q);
%             [bino_vector(outcome) outcome layer_length rank]
        end
    end
        
else
    Pr=0;
end

assert(isnan(Pr)==0,['MatrixFun returned NaN with:' num2str([pkts_recv,layer_length,rank,g,q])]);

end

% Should work
function PMWR = ProbMatricesWithRank(m,n,r,q)

% Get first set of gaussian coefficients
gc=gausscoeffs2(n,r,q);

% % % Calculate "sum"
% % val=0;
% % for k=0:r
% %     % This should be the one!
% %     val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(m*k+binomcoeffs(r-k,2)-n*m));
% % end

% Calculate "sum"
val=0;
for k=0:r
    % This should be the one!
    val1=gausscoeffs2(r,k,q);
    val2=q^(m*k+binomcoeffs(r-k,2)-n*m);
      
%     assert(isnan(val1)==0,'val1')
%     assert(isnan(val2)==0,'val2')
    
    % Multiply by 0 before inf :)
    if val1==0 || val2==0
        comb_val=0;
    else
        comb_val=val1*val2;
    end
    
    val=val+((-1)^(r-k)*comb_val);
    
%     assert(isnan(val)==0,'ProbMatricesWithRank returned NaN');
    
end

% Return probability of matrix 'm'x'n' with rank 'r'
PMWR=gc*val;

assert(isnan(PMWR)==0,['ProbMatricesWithRank returned NaN, params: ' num2str([m n r q])]);

end

% Should work (Tested! see bottom)
% Deprecated
function GC = gausscoeffs(m,r,q)
if r==0
    % disp('r = 0 in gauss coeffs')
    GC=1;
elseif r>0
    % disp('r > 0 in gauss coeffs')
    
    % Calculate numerator
    num=1;
    for w=m:-1:m-r+1
        num=num*(q^w-1);
    end
        
    % Calculate denominator
    denom=1;
    for w=r:-1:1
        denom=denom*(q^w-1);
    end
   
    % Calculate gaussian coefficient
    GC=num/denom;
    
elseif r<0
    disp('r < 0 error in gausscoeffs!!!')
end


end

% Should work (Tested! see bottom)
% Numerical improvements over gausscoeffs
function GC = gausscoeffs2(m,r,q)
if r==0
    % disp('r = 0 in gauss coeffs')
    GC=1;
elseif r>0
    % disp('r > 0 in gauss coeffs')
    
    % Calculate numerator
    % There is X factors in numerator
    
    
    num_length=abs((m))-abs((m-r+1));
    num=ones(num_length,1);
    num_index=1;
    
    for w=m:-1:m-r+1
        num(num_index,1)=(q^w-1);
        num_index=num_index+1;
    end
        
    % Calculate denominator
%     denom=1;
    
    denom_length=r;
    denom=ones(denom_length,1);
    denom_index=1;
    
    for w=r:-1:1
        denom(denom_index,1)=(q^w-1);
        denom_index=denom_index+1;
    end
       
    % Calculate gaussian coefficient
    GC=prod(num(:)./denom(:));
    
    assert(isnan(GC)==0,'Exception: gausscoeffs2 returned NaN')
    
elseif r<0
    disp('r < 0 error in gausscoeffs!!!')
end

end

% Not tested!
function bc = binomcoeffs(a,k)
% As on page 123 in "A course in combinatorics"

tmp_vector=ones(2,1);
tmp_index=1;

for w=0:-1:-k+1
    tmp_vector(tmp_index)=(a+w);
    tmp_index=tmp_index+1;
end

num=prod(tmp_vector);

denom=factorial(k);
bc=num/denom;

end




